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Monte Carlo simulations of the directional-ordering transition in the two-dimensional classical and 

quantum compass model 
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A comprehensive study of the two-dimensional (2D) compass model on the square lattice is performed for 
classical and quantum spin degrees of freedom using Monte Carlo and quantum Monte Carlo methods. We 
employ state-of-the-art implementations using Metropolis, stochastic series expansion and parallel tempering 
techniques to obtain the critical ordering temperatures and critical exponents. In a pre-investigation we re- 
consider the classical compass model where we study and contrast the finite-size scaling behavior of ordinary 
periodic boundary conditions against annealed boundary conditions. It is shown that periodic boundary con- 
ditions suffer from extreme finite-size effects which might be caused by closed loop excitations on the torus. 
These excitations also appear to have severe effects on the Binder parameter. On this footing we report on a 
systematic Monte Carlo study of the quantum compass model. Our numerical results are at odds with recent 
literature on the subject which we trace back to neglecting the strong finite-size effects on periodic lattices. The 
critical temperatures are obtained as Tc = 0.1464(2) J and Tc — 0.055(1)J for the classical and quantum 
versions, respectively, and our data support a transition in the 2D Ising universality class for both cases. 

PACS numbers: 02.70.Ss, 05.70.Fh, 75.10.Jm 



I. INTRODUCTION 

The compass model is one of the simplest models pos- 
sessing orbital degenerate states. Originally developed^ as 
a model for Mott insulators it has recently seen renewed 
interest*^ in connection with orbital-order in materials like 
transition metal (TM) compounds. Despite its closeness to 
ordinary models of quantum magnetism, like the Heisenberg 
model, there is no ordered phase characterized by magneti- 
zation properties. This means that the ordered phase appear- 
ing in the model is especially interesting in that it cannot be 
classified according to the Mermin- Wagner criterion. A com- 
petition of interactions in different directions rather results 
in a special long-range ordered state^ possessing a sense of 
orientation,-'' and the transition is at the same time accompa- 
nied by dimensional reduction.- The current interest in this 
model is furthermore triggered by the recent discovery that it 
describes arrays of superconducting Josephson junctions and 
because of a possible realization of a system which protects 
qubits against unwanted decay in quantum computation^-^ 

The compass model is a spin model on simple-cubic lattices 
in d dimensions of size N = L'^ defined by the Hamiltonian 



and the Hamiltonian assumes the form 



^ = (1/4) }2 (^-«+e. + J^^^^^+e. 



(2) 



^ — 2^ 2^ JkSi s'j_ 



+ek 



(1) 



where S^ represents the k-th component of a spin S at site 
i and i + e^ is the nearest neighbor of i in the k direction. 
In the classical case we have S G 0{d), or in a more ex- 
plicit vector representation with (p and 9 being angles on the 
sphere, we use the expression S^ = (cos((/3),sin((^)) and 
S^ = {cos{(p) sin(6'), sin((^) sin(6'), cos(6')) in two and three 
dimensions, respectively. In the two-dimensional (2D) quan- 
tum case S represents a spin-1/2 operator S = (1/2) (ax, Oz) 



where we have chosen the z instead of the y direction as a 
matter of convenience (usually we take S^ as the quantization 
component in quantum Monte Carlo). In this work the cou- 
pling constants are taken to be equal, Jk ~ J , and positive 
although the sign plays no role since it can be transformed 
away on bipartite lattices (L must be even). 

Recent contributions in the literature have explicitly inves- 
tigated the properties of the 2D compass model for both the 
classical and quantum Hamiltonian. Analytical and Monte 
Carlo work on the classical case proved the existence of 
a directional-ordering transition at finite-temperatures and it 
was argued that this transition belongs to the 2D Ising uni- 
versality classJ. Using exact diagonalization techniques and 
Green-function Monte Carlo the energy spectrum of low ly- 
ing states was analyzed for the quantum model in detail^i^ 
These studies provided the key result that the ground state is 
exponentially degenerate possessing a degeneracy of 2 x 2^. 
This turns the relatively simple Hamiltonian into a hard prob- 
lem comparable to frustrated magnets. Later work'- deter- 
mined the nature of the quantum phase transition to be of first 
order when driving the system by changing the coupling ra- 
tio Jxl Jz- A variant of the model possessing a similar quan- 
tum phase transition was finally analyzed in one dimension Jl 
In a recent Letter— the finite temperature properties of the 
quantum compass model were analyzed for the first time by 
means of a world line quantum Monte Carlo scheme based on 
the Suzuki-Trotter discretization. The authors conclude with 
the intriguing effect, that the presence of random site dilution 
has much weaker effects on criticality for quantum degrees of 
freedom than for classical ones. The numerical analysis sup- 
porting this conclusion is, however, based on rather small lat- 
tice sizes and the quality of the quantitative results is modest 



and in view of our results reported below would need further 

investigations. 

Due to the relevance of the model and the potential im- 
plications for future applications it would be desirable to 
have a more precise understanding of the critical behavior 
at the directional-ordering transition in the quantum compass 
model. The purpose of this work is to tackle this problem with 
a comprehensive Monte Carlo study for both the classical and 
quantum case where we will focus here on the non-disordered 
case. Our motivation to restudy the classical case is to gain 
as much experience as possible about the transition and diffi- 
culties that may arise in the Monte Carlo sampling and data 
analysis. Using this experience a large-scale simulation of the 
quantum compass model in 2D will follow in the second part. 
The next section introduces the methods and tools we used to 
accomplish this. Section|III]contains our results for the classi- 
cal compass model and Sec.|IV]the respective analysis for the 
quantum case. We close in Sec. |V] with a summary and our 
conclusions. 
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II. OBSERVABLES AND METHODS 



A. Observables 



In this section we describe the observables that are used 
to characterize the phases and to probe the phase boundaries 
of the compass model. The basic quantity is the total energy 
E = H = J2i, Ek and the corresponding heat capacity C = 
dE/dT. With Ek = JkS^S^j^^^ we denote the energy along 
the fc-th direction or on fc-bonds in the system. Using this 
definition a useful order parameter in 2D can then be defined 



D = — 

N 



Ef J QX QX T QZ QZ > 



(3) 



— \E^ 



E, 



where N ~ L'^. The quantity D measures the excess energy 
in one direction compared to the other direction. If D > 
the system is said to possess long-ranged orbital or directional 
order whereas for D ^ the system is disordered. An alter- 
native definition for the order parameter 



D' = — ( min Ek 



d=2 



T.E'^/d 



(4) 



fe=i 



can be used to give a visualization and characterization of the 
different phases as in Fig. [1] On the lattice we thereby mark 
all bonds which have less than the average bond energy (those 
that contribute most to the partition function) and look at the 
global structure of the resulting bond clusters. In the disor- 
dered phase we expect rather random clusters whereas the or- 
dered phase is characterized by clusters which are direction- 
ally ordered and independent of each other (dimensional re- 
duction). Note, that in two dimensions D and D' are actually 
the same quantity up to a constant factor, because D = 2D'. 
However, Eq. ^ provides the general possibility to define an 
order parameter in any dimension d, which might be useful 
for future studies. In order to investigate the universality class 
of the phase transition we further look at the susceptibility x 
and Binder parameter Q2 which are respectively defined as 



X = N{{D')-{Df), 



1.^1-^^'^ 



3(Z?2) 



2\2 



(5) 



where (D") denotes an average of the n-th moment computed 
from the time series of D. 

For the susceptibility we expect a finite-size scaling behav- 
ior of the form 
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FIG. 1: (Color online) Visualization of different phases in the 2D 
compass model. Left: For T > Tc the system is disordered and the 
distribution of bonds possessing less than average bond energy (thick 
lines) is rather random. Right: For T < Tc the prevalent correlations 
order into one direction, i.e. the system is in a directionally-ordered 
state. The pictures are snapshots of a Monte Carlo simulation of the 
classical model with L = 12 at T = 0.3 J and T = 0.10 J respec- 
tively (ferromagnetic representation). The small arrows indicate the 
spin degree of freedom. 



at the critical point with v being the correlation length critical 
exponent and 7 the exponent for the susceptibility. Neglect- 
ing corrections to scaling, the Binder parameters for different 
lattices sizes L should ideally cross at the critical temperature 
Tc- In any case, the behavior of the crossing points of Q2 at 
lattice sizes L and 2L should approach Tc like L^^''^^^ if we 
have corrections to scaling (w < cx)). 



B. Monte Carlo methods 

Ordinary Metropolis Monte Carlo simulations are used for 
the classical model where we update each spin sequentially. 
During the thermalization procedure we adjust the proposed 
moves such that an average acceptance rate of about 50% is 
obtained at each temperature. As it already becomes apparent 



from simulations on very small lattice sizes L that the system 
suffers from huge autocorrelation times we add a parallel tem- 
pering (PT) schemeliii, where we propose to exchange spin 
configurations between simulation threads at different temper- 
atures Ti. This exchange is attempted every n sweeps, where 
n is typically in the range 2 to 20. By tracking individual con- 
figurations we make sure all temperatures are seen and that 
sufficient diffusion through temperature space is performed. 
For simplicity the simplest PT scheme is used, meaning that 
an equidistant temperature spacing between neighboring pro- 
cesses is chosen. As a result a reduction of autocorrelation 
times by two orders of magnitude is achieved which pays off 
in comparison to little longer simulation times. 

In case of quantum spin degrees of freedom, we employ a 
quantum Monte Carlo (QMC) procedure based on the stochas- 
tic series expansion (SSE)JS technique originally developed 
by Sandvik. Our own implementation is based on the (di- 
rected) loop schemei- supplemented by ideas of Ref. [O- Re- 
call that the principle of SSE is sampling the series expansion 
of the quantum partition function 



z^tr{cM~m))^Y.Y. 
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by a Markov chain stochastic process, where /3 ~ l/fceT is 
the inverse temperature. The last line of Eq. ^ is the central 
starting pointi^ of the method because it specifies the config- 
uration space (and the weights) in which the sampling takes 
place. A configuration lives in the product space of spin con- 
figurations \a) times the space of all possible sequences (or 
permutations) Sn of n bond operators (or vertices) Hbf ■ The 
degrees of freedom are thus \a), n, and Sn, which are sam- 
pled by the usual combination of diagonal, non-diagonal, and 
spin flip updates."' 

In the case of the compass model the bond operators Ti),. 
can be derived from the Hamiltonian ^ as 



allowed. Unfortunately, this asymmetry in the operator rep- 
resentation cannot be transformed away by a simple "sym- 
metrizing" rotation of the Hamiltonian because of emerging 
minus sign problems. Note finally that the non-zero energy 
shift e has an effect on the order parameter D since it influ- 
ences the number and the distribution of bond operators in the 
operator sequenced This effect can cause additional finite- 
size contributions also in the susceptibility and the Binder 
parameter which vanish in the thermodynamic limit and for 
T — *■ 0. We have checked that at L = 16 no difference 
could be detected in the susceptibility maxima locations for 
e = 0.1, 0.5, 0.9 within error bars. Here we work with e = 0.5 
at all lattices sizes. 

Since simulations of the quantum model display the same 
rapid critical slowing down as the classical model we perform 
additional quantum PT updates^S, in the same manner as de- 
scribed above. We implemented both the classical and quan- 
tum Monte Carlo PT scheme on parallel architectures with 
the restriction of shared memory access for fast communica- 
tion between processes. This is essential since PT updates are 
done rather often. 

For data analysis purposes we use well-known multi- 
histogram techniques to optimally combine simulations at 
different temperatures. Those techniques are available for 
both the classical'^ and quantum cases22iSi. In combination 
with optimization routines like the Brent method^^ they al- 
low rather systematic and unbiased estimation of pseudocriti- 
cal temperatures from peaks in the susceptibility. 



C. Boundary conditions 

Ordinarily, the vast majority of Monte Carlo simulations are 
performed using periodic boundary conditions (pbc) which 
map the lattice onto a torus topology using the assumption 
that free-energy contributions from the surface are thereby 
minimized. In contrast to this approach, Mishra et al^ ar- 



HbG 



I SfS', 



if 6 is a z bond 



I { S+ S+ , Sr SJ , S+ SJ- ,SrS+} if 5 is a X bond , 



where the appearance of pure S^ S* , and S^ S^ terms are a 
notable difference to an ordinary Heisenberg model. Here S^ 
and S~ refer to creation and annihilation operators and the 
subscripts i,j are the two sites of the bond b. Simulations of 
the quantum compass model are furthermore more involved 
since the Hamiltonian dictates an asymmetry between bonds 
in X and z direction, allowing no spin flip operators of type 
S^S^ to reside on z-bonds. On the other hand, there are a 
priori no diagonal terms S^S^ on .x -bonds. However, since 
non-diagonal terms can only be introduced into the SSE con- 
figuration space after the diagonal-update (non-diagonal op- 
erators must be present!) we are therefore forced to introduce 
a positive non-zero energy shift e into the Hamiltonian of Eq. 
(|2|i. As a consequence both non-diagonal and diagonal terms 
may reside on .i-bonds. On z-bonds only diagonal terms are 



pbc 



abc 



FIG. 2: (Color online) Visualization of the different boundary condi- 
tions used in this work. Left: Ordinary periodic boundary conditions. 
All bonds carry the same coupling and the dashed bonds connect 
the spins across boundaries. The topology is a torus. We refer to 
this case as pbc. Right: So called "annealed" boundary conditions 
(abc).^ Here the sign of the couplings on the dashed boundary bonds 
may fluctuate dynamically resulting in an additional degree of free- 
dom. As an example we draw some thick bonds indicating a negative 
coupling. 
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FIG. 3: Data for the 2D classical compass model obtained by Monte Carlo simulations. The top row displays the results for periodic boundary 
conditions and the bottom row for annealed boundary conditions. Note the the temperature ranges are different for both cases and that not all 
lattice sizes are shown for better readability. The lines through the data points are obtained from the multi-histogram analysis. Left: The order 
parameter D as a function of different lattice sizes L. Middle: The susceptibility x of the order parameter. Right: The Binder parameter Q2 . 



gue in their recent contribution that periodic boundary condi- 
tions might not be optimal in the case of the compass model. 
Instead, they introduce special, so called annealed boundary 
condition (abc) to arrive at their Monte Carlo results. Since a 
detailed comparison between these two boundary conditions 
has to our knowledge not been done, we will explicitly study 
and compare their effect on the finite-size scaling behavior for 
the classical compass model. This comparison is especially 
interesting in view of the fact that we may not easily apply 
the annealed case to quantum Monte Carlo since it induces a 
minus sign-problem. A characterization and understanding of 
the scaling behavior for periodic boundary conditions would 
therefore be of advantage before studying the quantum case. 

Figure|2]displays these two types of boundary conditions as 
a sketch. The topology of the annealed boundary condition is 
the same as for periodic boundary conditions. The annealed 
case is special because the sign of couplings on bonds across 
the border may fluctuate dynamically according to the Boltz- 
mann distribution. The bond sign is therefore an additional 
degree of freedom in the Monte Carlo update rendering the 
simulations somewhat more complex. 



III. THE CLASSICAL COMPASS MODEL IN 2D 

In this section we start the presentation of our simulation 
results. We consider firstly Monte Carlo simulations of the 
2D classical compass model. The main purpose of this sec- 
tion is to give an explicit comparison between the different 
boundary conditions introduced in the last section. To this 
end we run simulations for both cases and compare the ob- 
servables of Sec. |ll] and their finite-size behavior. Figure [5] 
gives an overview of our Monte Carlo estimates for D, x and 



Q2- There, the top row contains results for periodic boundary 
conditions and the bottom row for annealed boundary condi- 
tions. Knowing the different behavior of reaching the thermo- 
dynamic limit is useful in order to appreciate results of our 
simulations which follow in subsequent sections. 

Simulations are done using lattice sizes L = 
{10,12,16,24,32,36,48,64,128} periodic be and 
L = {10, 12,16,20,24,36,52,64} for annealed be, typ- 
ically taking about 10^ measurements per data point after 
an equilibration phase of 10^ sweeps. By the behavior of 
the order parameter in Fig. [5] (left) it is immediately evident 
that there is a phase transition and that directional order with 
Z) > is realized at low temperatures. We secondly observe 
that the order parameter for the pure periodic case has a slow 
convergence for small lattice sizes while for larger sizes it 
suddenly moves considerably. In contrast, the data for the 
annealed case show a much smoother movement towards the 
infinite-volume limit and it is evident that finite-size effects 
are drastically reduced. A difference like this is actually 
expected for different boundary conditions. The crucial and 
interesting question is whether the two boundary conditions 
lead to the same critical temperature in the infinite volume 
limit where boundary effects should vanish. 

We therefore obtain an estimate of the critical point T^ in 
the thermodynamic limit by fitting the pseudocritical temper- 
atures Tc{L) taken from the peaks of the susceptibilities in 
Fig.[3](middle) at lattice size L to the finite-size scaling ansatz 



Tc(i) = Tc + bL-'^/"{l + cL 



(8) 



Here b, c are some constants and u; is an exponent describ- 
ing corrections to scaling. In a first step, we assume nothing 
about the value for the correlation length exponent ly and leave 
it as fit parameter. The fitting procedure to the data in Fig. |4] 
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FIG. 4: (Color online) Determination of the critical temperature 
Tc from finite-size scaling of the pseudocritical temperatures de- 
termined from susceptibility peaks. The two curves correspond to 
different boundary conditions which trigger a completely different 
convergence to the critical point. The lower curve is obtained for 
annealed boundary conditions (abc) and shows the superior scaling 
compared to periodic boundary conditions (pbc). Lines are fits to 
Eq. ([8} neglecting the correction term uj. 



yields Tc = 0.144(2) J from periodic boundary conditions 
and Tc = 0.1461(8) J from annealed boundary conditions 
and the estimate for the correlation length critical exponent 
is I' = 0.98(4) which we take from the straight line fit for 
the annealed case. Both results agree within error bars. The 
annealed value yields a much more accurate estimate since 
here the asymptotic scaling regime sets in much earlier and 
we have more points available for fitting. These numerical 
estimates are in accordance with the value Tc = 0.147(1) J 
obtained in Ref. S With our value for v we support the claim 
that the transition is of 2D Ising type. To further confirm this 
conjecture we also determine the exponent 7 associated with 
the susceptibility x- For lattice sizes large enough {L > 20) 
we obtain j/iy = 1.73(4) from the annealed case (see Ta- 
bleUand Fig.|2]below) which is again consistent with 2D Ising 
universality. In a second step, we can now assume Ising uni- 
versality to be given to improve the fit. Using ly = 1 as a 
fixed parameter the improved value for critical temperature is 
Tc = 0.1464(2)J. 

Let us now turn to a discussion of the Binder parameter 
Q2 displayed in Fig. [5] (right). For the annealed case a nice 
crossing of curves at the critical temperature can be observed 
and our estimate for the Binder parameter at the crossing 
point (taking the three largest lattice sizes) is Q2 = 0.61(1). 
This is roughly the known value for the 2D Ising model, 
which - however - is usually obtained for periodic bound- 
ary conditionsi^ Using the observed crossing behavior, the 
Binder parameter supplies a natural third check of the crit- 
ical temperature and the critical exponent. We hence ap- 
ply our recently developed data collapsing tool-- and obtain 
Tc = 0.1465(4) J and 1/ = 1.01(4) from the best data col- 
lapse. These values are again fully consistent with our re- 
sults above and give further confidence to our analysis. In 
contrast, the nice properties of the Binder parameter do not 



show up for periodic boundary conditions, where it is hard to 
judge whether curves for different lattice sizes cross in a sin- 
gle point at all. Rather, we see strong finite-size effects and 
that the crossing points for large lattice sizes move close to 
2/3, which is totally in contrast to the expected behavior It is 
known that different boundary conditions cause a discrepancy 
(see for instance Refs.l23landl25h in the Binder crossings but 
such a drastic behavior was unexpected. 

In summary, our investigation for the classical model 
clearly show that annealed boundary conditions are favorable 
because they drastically reduce finite-size effects and yield 
good scaling properties for the finite-size analysis. With pe- 
riodic boundary conditions much larger lattice sizes need to 
be investigated in order to obtain the critical temperature and 
to approach the right asymptotic scaling regime. Addition- 
ally, our analysis shows that the Binder parameter for periodic 
boundary conditions does not cross at the usual expected value 
and that we may not use the crossing point (height) as a good 
indication for the critical point, whereas for annealed condi- 
tions we get good properties. These effects are currently not 
properly understood. By referring to the typical spin configu- 
ration in Fig.[T]it is, however, tempting to argue that the dom- 
inant energy correlations (blue lines) wrap around the torus in 
the ordered phase thereby forming some kind of closed loop 
excitations. These excitation appear to be more stable against 
thermal fluctuations than open excitation. Annealed bound- 
ary conditions seem to prohibit the formation of such loops 
leading to a better scaling behavior 



IV. THE QUANTUM COMPASS MODEL IN 2D 

Using the knowledge gained from simulations of the clas- 
sical compass model we turn to the discussion of the simula- 
tion results of the quantum version. Simulations are done us- 
ing the stochastic series expansion as outlined in Sec. HIl The 
reader is reminded that annealed boundary conditions, where 
the sign on boundary bonds fluctuates, are not possible be- 
cause such fluctuations induce a sign problem in the quantum 
Monte Carlo scheme. 

We therefore choose to simulate with periodic bound- 
ary conditions and expect from Fig. |4] that large lat- 
tice sizes might be needed to see the right scaling 
and to obtain the infinite-volume critical temperature. 
Using the parallel tempering scheme and the reduc- 
tion of autocorrelation times by two orders of mag- 
nitude, we were finally able to simulate lattice sizes 
L = {8, 10, 12, 14, 16, 18, 20, 24, 28, 32, 40, 48, 52, 64}, 
where the largest one is about the limit one can reach in 
quantum Monte Carlo in feasible time and resources at the 
momenti^^ Our largest system size is about three times as 
large compared to the simulations of Ref. [Tj- A detailed 
check and verification of our algorithm was done with data 
from full exact diagonalization (ED) on a 4 x 4 lattice. We use 
our own ED program with some implemented symmetries^ 
to reduce the dimension of the Hilbert space, as well as the 
ALPS package^ for smaller system sizes. During the Monte 
Carlo runs, a total number of about 4 x 10^ measurements 
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FIG. 5: QMC results for the 2D quantum compass model 
with periodic boundary conditions. All lines are a guide 
to the eye. (a) The order parameter D for lattice sizes 
L = {12, 14, 16, 18, 20, 28, 32, 40, 52, 64} displays a clear signal of 
a stable ordered phase at low ternperatures. The arrow marks the tran- 
sition temperature Tc from Ref.[l2|. Our own data indicates a smaller 
value, (b) The susceptibility x on a logarithmic scale for lattice sizes 
L = {10, 12, 14, 16, 18, 20, 28, 32, 40, 52, 64}. (c) the Binder pa- 
rameter Q2 in the quantum compass model with periodic boundary 
conditions, where steeper slope corresponds to larger L (neglecting 
L = 28 and L = 48 for better clarity). The qualitative behavior is 
the same as for the Binder parameter with periodic boundary condi- 
tions in the classical model. No common crossing point is present 
for the lattice sizes of this work. 
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FIG. 6: (Color online) Finite-size scaling of pseudocritical tempera- 
tures for different lattice sizes obtained from the susceptibility. For 
small lattice sizes corrections to scaling are evident. For large lattice 
sizes 2D Ising scaling is reached yielding our estimate for the criti- 
cal temperature of Tc/J — 0.055(1) J. The curve trough the points 
represents a fit to Eq. ^. We also show the crossing points of the 
Binder parameters at L and 2L for a consistency check. The arrow 
indicates the previous result of Ref. 12. 



are typically taken after each sweep and 2 x IC* sweeps are 
used for thermalization. Those numbers are, of course, only 
meaningful with the additional information that we construct 
as many loops in the non-diagonal update such that on aver- 
age 2n vertices are visited in the SSE configuration. Figure|5] 
shows the result for the order parameter, the susceptibility and 
the Binder parameter obtained from the simulations in this 
manner The behavior of the order parameter shows a clear 
signal of a transition from a disordered to an ordered state 
at small temperature, evidently becoming more pronounced 
with increasing lattice size. This proves the existence of a 
directional-ordering transition also in the quantum case. In 
the ordered phase, the order parameter seems to take on a 
value which is quite different from the classical case and the 
order is furthermore less stable against thermal fluctuations 
as the temperature regime in the quantum case is evidently 
much smaller Note that the overall estimate for D also agrees 
roughly with data of Ref.[l2l 

The dependence of the data on the lattice sizes is, as ex- 
pected, qualitatively similar to the classical case, i.e., the or- 
der parameter curves and the susceptibility peaks shift con- 
siderably to lower temperatures for larger lattice sizes. This 
shift is in fact so large that it is already obvious from Fig. |5] 
(a) that the previous estimate of the critical temperature in the 
literature is much to large. Before we quantify this discrep- 
ancy for the critical temperature, we draw our attention to the 
susceptibility and the Binder parameter in Figs.|2b),(c), both 
showing a behavior similar to the classical case with periodic 
be. We note especially that the Binder parameter is again be- 
having oddly and that there is not a well defined crossing seen 
at all at the lattice sizes simulated. A crossing point might 
still be achieved for very large lattice lengths L but is cer- 
tainly difficult to quantify since the value of Q at the crossing 
point is very close to 2/3. Due to this observation the Binder 
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FIG. 7: (Color online) Plot of the susceptibility x (tit peak locations) 
versus system size L for all different simulations in this work (clas- 
sical annealead be, classical periodic be, quantum periodic be in this 
order from top to bottom) on a double logarithmic axis. The straight 
lines are fits to Eq. l|6). All cases are consistent with a value of 
"ilv = 1.75. 



parameter is clearly not suited to determine the critical tem- 
peratures by looking at the Binder crossings for small lattice 
sizes, where the true behavior is just not seen. 

Let us now determine an improved value for the critical 
temperature with finite-size scaling from the maxima of the 
susceptibilities. As in the classical case, we fit the pseud- 
ocritical values to the scaling ansatz given in Eq. ^. To 
use as many data points as possible we include coiTections to 
scaling, described by the exponent cj, into the fit and leave 
all fit parameters free. Including all lattice sizes we ob- 
tain Tc = 0.055(6) J and v = 0.9(2) with a fit quaUty of 
X^/d.o.f = 0.66. Those values, however, are stable also for 
fitting windows starting at larger lattice sizes. The precision 
for the critical exponent v is rather low but agrees with our 
expectation of 2D Ising universality within the error bar. Un- 
der this assumption, we fix j/ = 1 and repeat the fit procedure 
yielding an improved estimate for the critical temperature as 
Tc = 0.055(1) J. The relative discrepancy with the previous 
estimate of Ref. |l2J is approximately 36%. As a cross check 
for our analysis we further look at the scaling of the crossing 
points of Q2 at lattice sizes L and 2L, which is also indicated 
in Fig. |6l We observe that this scaling is consistent with the 
previous value from the susceptibilities but we do not attempt 
a detailed fit by lack of enough data points. It is then also 
useful to obtain an independent estimate of the critical tem- 
perature from the maxima in the heat capacity C, which again 
gives consistent results but does not reach the accuracy of our 
previous analysis since C is generally hard to sample in QMC 
at low temperatures. 

To finalize our analysis, we determine the critical exponent 
7 from the susceptibility of the order parameter D. For large 
lattice sizes we expect a scaling according to Eq.|6] which can 
be tested by plotting ln(x) versus ln(L). This is done in Fig.]?] 
together with the data for the classical cases. It is evident that 
asymptotic scaling sets in only for the largest lattice sizes from 
which we obtain a value of 7/1/ — 1.68(8) consistent with 



TABLE 1: Results for the critical temperature and critical exponents 
as obtained in this work. The upper section contains the results for 
the classical model taken from annealed be, which are all pairwise 
consistent. The middle section summarizes our estimates for the 
quantum compass model for the cases with and without the assump- 
tion of 2D Ising behaviour {v — 1). Both cases are consistent with 
each other. Lastly the values for ^/v are summarized as obtained 
from the largest lattice sizes for the different simulation runs. 
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LO 
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no assumption 
2D Ising 
collapse Q2 


0.1461(8) 
0.1464(2) 
0.1465(4) 


0.98(4) 
1.01(4) 


- 


1.3 
1.15 


no assumption 
2D Ising 


0.055(6) 
0.055(1) 


0.9(2) 


0.5(4) 
0.8(2) 


0.66 
0.61 




class, (abc) 


class, (pbc) quant, (pbc) 




j/iy 


1.73(4) 


1.72(5) 


1.68(8) 





2D Ising universality, but not precise enough to be absolutely 
conclusive. 



V. SUMMARY AND CONCLUSIONS 

In this paper we reported on comprehensive Monte Carlo 
simulations of the classical and quantum compass model. 
By comparing different boundary conditions for the classical 
case, we showed that for ordinary periodic boundary condi- 
tions one needs to go to very large lattice sizes to see the right 
scaling and to get good convergence to the critical point. In or- 
der to simulate large lattice sizes, we implemented a parallel 
tempering scheme to counteract huge autocorrelation times. 
Our results, which are summarized in Table H] are perfectly 
consistent with previous studies in the literature for the classi- 
cal model. For the quantum model our simulations are quanti- 
tatively at odds with earlier studies and we provide here a new 
estimate for the critical temperature Tc- We argued that this 
discrepancy might be explained by the huge finite-size cor- 
rections originating from stable loop excitations formed by 
correlation orderings which appear on the torus topology at 
periodic boundary conditions. It appears that those excita- 
tions even destroy the usual properties of Binder parameters. 
Our analysis, however, shows that one can still arrive at an 
estimate for Tc at periodic boundary conditions provided that 
one takes this effect into account. All critical exponents ob- 
tained in this study give further support to the claim that 2D 
Ising universality describes the directional-ordering transition 
in the 2D compass model. 

Our findings for the quantum model might have an impact 
on the conclusions of Ref. [Ij because a precise estimate of 
Tc enters into the analysis of dilution effects on the model. 
It is conceivable that the conclusion obtained there are still 
qualitatively valid. For a true quantification of the dilution 
effect, however, there is no way around performing a more 
detailed investigation of larger lattice sizes. The knowledge 
gained in this work should help to start such a study. 
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The precision of our results for the quantum model are 
still rather low compared to many other systems of statisti- 
cal physics. In this respect it would be an interesting future 
project to devise and analyze special boundary conditions for 
the quantum model with improved finite-size scaling behavior 
compared to periodic boundary conditions. 
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